module emit

  use globals
  use space
  implicit none
  
contains
!--------------------------------------------------------------------------

subroutine read_qnf

  implicit none
  integer  :: ii, jj, qnx, qny
  
  open(25,file=trim(adjustl(ipath))//trim(adjustl(nf_in)), form='FORMATTED',status='old')
!  qny direction = azimuthal (should be 2pi total)
!  qnx direction = polar  (should be pi total)
  read(25,*) qnx, qny
  
!  write(*,*) "Size of Qnf field: ", qnx, qny
  allocate(Qnf(qnx,qny))

!  if: rows = azimuthal
!      columns = polar 
!  read(25,*) ((Qnf(ii,jj),ii=1,qnx),jj=1,qny)

!  if: rows = polar
!      columns = azimuthal 
  read(25,*) ((Qnf(ii,jj),jj=1,qny),ii=1,qnx) ! effectively transposing

  close(25)
  
  Qtot = trap3(Qnf,qnx,qny) ! Form: Qtot = trap3(Qnf(p,a),polar,azimuthal)
!  print *, "Qtot: ", Qtot

end subroutine read_qnf
!--------------------------------------------------------------------------

subroutine read_fowler

  implicit none 
  integer  :: ii
  
  open(23,file=trim(adjustl(ipath))//trim(adjustl(fo_in)), form='FORMATTED',status='old')
  read(23,*) nchi
  
  allocate(rchi(nchi))
  allocate(rf(nchi))
  
  do ii=1,nchi
     read(23,*) rchi(ii), rf(ii)
  enddo
  
  close(23)

end subroutine read_fowler
!--------------------------------------------------------------------------

function fowler(chi)

  implicit none
  real(8) :: fowler, chi
  integer :: nlo, nhi, ii
  
  
  if (chi .lt. rchi(1)) then
     fowler = 0.0d0
  elseif (chi .gt. rchi(nchi)) then ! extrapolating here, not the best...
     fowler = (rf(nchi)-rf(nchi-1))/(rchi(nchi)-rchi(nchi-1)) * &
	        & (chi-rchi(nchi-1)) + rf(nchi-1)
  else ! linear interpolation
     do ii=1,nchi
	    if (chi .gt. rchi(ii)) then
	       nlo = ii
		   nhi = ii+1
		   exit
		endif
	 enddo
     fowler = (rf(nhi)-rf(nlo))/(rchi(nhi)-rchi(nlo)) * &
	        & (chi-rchi(nlo)) + rf(nlo)
  endif
  
end function fowler
!--------------------------------------------------------------------------

!  Too much to explain here, see my thesis...
!  Note: this is HARD-CODED FOR 3-photon process
subroutine emission

  implicit none
  real(8)  :: Eone
  
  Ao = 120.0d0*100.0d0**2.0d0
  CC = 1.0d0*10.0d0**(-7.0d0)/((1.0d0*10.0d0**6.0d0)**3.0d0 * 100.0d0**4.0d0)
  
  ! note: mistake in thesis here, forgot to include As in c_therm
  c_therm = 4.0d0*pi*rad**3.0d0*Ao
  c_mpe = 2.0d0*CC*rad**3.0d0*(2.0d0*Fl/tau*(log(2.0d0)/pi)**(1.0d0/2.0d0))**3.0d0
  
  chi = (ord*h*(c/w)-bh*eVJ)/(kb*Te)
  Fo = fowler(chi)
  
  ftherm = Te**2.0D0 * exp((-bh*eVJ)/(kb*Te))
  fmpe = Fo/chi**2.0D0 * exp(-12.0D0*log(2.0D0)*((timep-2.0D0*tau-toffset)/tau)**2.0D0)
  
  therm = therm + dtp*(ftherm+ftherm_old)/2.0D0
  mpe = mpe + dtp*(fmpe+fmpe_old)/2.0D0
  
  Efield = 1.0D0/(3.0D0*Vp*eo) * (c_therm*therm+c_mpe*mpe*Qtot)  ! induced efield at particle surface, continuous
  Eone = 1.0d0*ec/(3.0D0*eo*Vp/rad)                              ! Efield to emit single electron
  etot = floor((Efield * 3.0D0*eo*Vp/rad)/ec)                    ! emitted electrons
  Efield = etot * Eone                                           ! corrected efield for discrete emission

  if (Efield .ge. Eone .and. eflag3 .eq. 0) then
     write(*,'(a,d11.5,a)') " >>> Electron emission begun at time ",timep," <<<"
	 eflag3 = 1.0d0
  endif 
  
  if (Efield .ge. Ecrit .and. eflag1 .eq. 0) then
     write(*,'(a,d11.5,a)') " >>> Particle ablation begun at time ",timep," <<<"
	 eflag1 = 1.0d0
  endif  
    
  if (Efield .ge. Ecrit*rad/latspace .and. eflag2 .eq. 0) then
     write(*,'(a,d11.5,a)') " >>> Particle completely ablated at time ",timep," <<<"
	 eflag2 = 1.0d0
	 ! etot = etot_old  ! capping etot as the analysis is really silly past this point
     ! this may create problems with interpolating for plots and needs to be moved anyway
  endif
  
  detot = etot - etot_old
  
  ftherm_old = ftherm
  fmpe_old = fmpe
  etot_old = etot
  
end subroutine emission
!--------------------------------------------------------------------------

function aeout_calc(dummy)

  implicit none
  
  real(8) :: aeout_calc
  real(8) :: eout, rdetot, aeout, aeoutd, aeoutn, ee
  real(8) :: aeoutd_old, aeoutn_old, ri, REN, maxef
  integer :: ii, EN, dummy
  
  EN = 100
  REN = EN
  maxef = 10.0d0
  aeoutd = 0.0d0
  aeoutd_old = 0.0d0
  aeoutn = 0.0d0
  aeoutn_old = 0.0d0

	  ! Esentially a center of mass integral over fermi-dirac distribution past barrier height  
	  ! a. energy deposited into water via electron emission
      !    this should be average energy of escaped electron less bh (done after function call)
      !    as calcualted from fermi-dirac distribution from bh to inf
	  ! b. energy removed from particle absorption via electron emission
	  !    this should be average energy of escaped electron including bh	  

	  do ii=1,EN+1
	     ri = ii
	     ee = bh*eVJ * (1.0d0 + maxef/REN*(ri-1.0d0)) ! THIS MAKES NO SENSE< DOUBLE CHECK!!!
	     aeoutd = aeoutd + (bh*eVJ*maxef/REN) * (aeoutd_old + 1.0d0/(exp((ee-Efermi*eVJ)/(kb*Te))+1.0d0))/2.0d0
	     aeoutn = aeoutn + (bh*eVJ*maxef/REN) * (aeoutn_old + ee/(exp((ee-Efermi*eVJ)/(kb*Te))+1.0d0))/2.0d0
	     aeoutd_old = aeoutd
	     aeoutn_old = aeoutn
      enddo
	  
	  aeout_calc = aeoutn/aeoutd

end function aeout_calc
!--------------------------------------------------------------------------

function ewater_heat(ii, we, weo, aeout)

  implicit none
  
  real(8)  :: ewater_heat, we, weo, dwe, mfp
  real(8)  :: nwater, aeout
  integer  :: ii
  
  nwater = rho/(molw*1000.0d0)  ! #/m^3
  dwe = we - weo
  mfp = 100.0d-9
  ewater_heat = (aeout*dwe)/mfp * exp(-1.0d0*r(ii)/mfp)

end function ewater_heat

!--------------------------------------------------------------------------

subroutine deallocate_emission

  implicit none
  
  deallocate(rchi)
  deallocate(rf)
  deallocate(Qnf)
  
end subroutine deallocate_emission
!--------------------------------------------------------------------------

end module emit
